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Abstract 

Nose and Hoover’s 1984 work showed that although Nose and Nose-Hoover dynamics were both 
consistent with Gibbs’ canonical distribution neither dynamics, when applied to the harmonic oscil¬ 
lator, provided Gibbs’ Gaussian distribution. Further investigations indicated that two independent 
thermostat variables are necessary, and often sufficient, to generate Gibbs’ canonical distribution 
for an oscillator. Three successful time-reversible and deterministic sets of two-thermostat motion 
equations were developed in the 1990s. We analyze one of them here. It was developed by Martyna, 
Klein, and Tuckerman in 1992. Its ergodicity was called into question by Patra and Bhattacharya 
in 2014. This question became the subject of the 2014 Snook Prize. Here we summarize the previ¬ 
ous work on this problem and elucidate new details of the chaotic dynamics in the neighborhood 
of the two fixed points. We apply six separate tests for ergodicity and conclude that the MKT 
equations are fully compatible with all of them, in consonance with our recent work with Clint 
Sprott and Puneet Patra. 
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I. DETERMINISTIC TIME-REVERSIBLE THERMOSTATS AND ERGODICITY 


In 1984 Shuichi Nose discovered a canonical dynamicsi*^ consistent with Willard Gibbs’ 
canonical phase-space distribntion. Hoover- used a generalization of Liouville’s flow equation 
to develop a “NosGHoover dynamics”, a simpler variation on Nose’s work. He pointed out 
that neither approach gave Gibbs’ complete canonical distribution for the simple harmonic 
oscillator problem : 

‘H{q,p) = [ {q^/2) + (pV 2) ] —t f{q,p) oc [ Gibbs ] . 

In the 1990s three generalizations of this work were developed^'- to remedy the stiffness 
and the lack of ergodicity that resulted when Nose’s ideas were applied to the harmonic 
oscillator. All three of them include two “thermostat” variables, { C)0 > which control the 
motion of the oscillator coordinate and momentum { q,p} , steering it toward the canonical 
distribution. The generalized motion equations all satisfy an analog of the hydrodynamic 
continuity equation, (dp/dt) = —V • (pu) . The stationary ( steady-state ) version of this 
phase-space flow equation is : 

(df/dt) = 0 = -{dfq/dq) - (dfp/dp) - {dfC/dQ - (dfi/d^ . 

For all three flow models the probability density f{q,PiCi€} is stationary when it includes 
Gaussian distributions for the two new thermostat variables, ( and ^ . 

As a result there are at present three sets of four ordinary differential motion equations 
all of which provide the full canonical distribution for an oscillator along with Gaussian 
distributions for the additional thermostat variables {C)0- 

{q = p ■ P = -q-(p-^p^ ■ ( = p^ -I ■ ^=p^ -Sp"^ } 

{q = p; p = -q - (^p - ^p^ ; ( = p"^ - I ■ ^ } [ JB ]^ ; 

{q = p- p=-q-(p- (=p^-l-^(- ^ = (^-1} [MKT]^. 

Each of them displays a “mirror” or “inversion” or “rotational” symmetry in the {q,p) plane: 
any solution { +q{t), +p(t), C{p)y^{'t) } has a mirror image when the oscillator is viewed in a 
mirror perpendicular to the q axis. The solution viewed in the mirror replaces both +q and 
+p by their mirror images, —q and —p . In the mirror solution, { —q(t), —p(t),({p),^(t) } > 
the time-dependent thermostat variables ( and ^ are unchanged. 
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There are also generalizations of each of these ideas based on controlling more, or dif¬ 
ferent, moments of the canonical distribution function.-i^ In addition, a variety of different 
solutions result for thermostat relaxation times other than unity, for coordinate-dependent 
temperature prohles, and for more complicated potentials. 

In addition to the canonical oscillator probability oc the thermostat variables 

(C,0 ^1^0 have Gaussian distributions : 

Jhh = Jmkt D ; fjB D e“‘’ . 

Patra and Bhattacharya^ investigated the {q,p) phase-space density in the vicinity of 
an unstable hxed point {q,p,CyO = (0, 0, —1,-|-1) of the MKT equations. They displayed 
an apparent low-probability region there and suggested that the MKT equations were not 
ergodic. Because any lack of ergodicity would contradict Martyna, Klein, and Tuckerman’s 
belief in the ergodicity of their own model, we established the 2014 Snook Prize^^ as a reward 
for the most convincing work demonstrating either ergodicity or its lack. In January 2015 
we awarded the prize to the authors of Reference 11. 

Here we clarify the differing conclusions of References 9 and 11 by exploring six aspects 
of the chaotic dynamics and stationary measure of the MKT equations. These include [A] 
the moments of the measure, [B] the largest Lyapunov exponent, [C] the two fixed points 
of the flow, [D] the attractor/repellor dynamics near the two fixed points, [E] the measure 
in the neighborhood of these hxed points, and [F] the symmetry of the measure in the 
neighborhoods of 81 lattice points arranged in a four-dimensional phase-space lattice. Our 
description of the underlying analysis and numerical work makes up the following Section 
II. Our Conclusions follow in the Summary Section III. 

II. INVESTIGATING ERGODICITY: THE MKT HARMONIC OSCILLATOR 

Ergodicity, with any dynamical trajectory coming close to all phase-space states, became 
an issue with the study of the one-thermostat Nose-Hoover oscillator-J^ : 

q = p ■ p = -q - (p ■ ( = p‘^ -I . 

This model exhibits a variety of regular solutions. Most trajectories correspond to two- 
dimensional tori in the three-dimensional {q,p,C) phase-space. About hve percent of the 
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Gaussian phase-space measure , 


makes up a chaotic sea perforated by the tori^i. 

Surprisingly, adding a fourth variable to the phase space has a tendency to simplify the 
flow, with the chaotic region expanding to £11 the entire phase space. In what follows we 
consider the details of the Martyna-Klein-Tuckerman oscillator : 

{q = p- p=-q-(p- ( = i = e-l} [MKT}^. 

All of the numerical work described here was carried out with the classic fourth-order Runge- 
Kutta integrator, mostly with a timestep of 0.001 . We consider six different aspects of the 
MKT oscillator’s phase-space flow, and show that all of them are fully consistent with the 
ergodicity of that model. 

A. Moments of the Distribution Function 

If MKT dynamics is ergodic then its long-time-averaged distribution is Gaussian : 

The independence of the four variables implies that the second, fourth, and sixth moments 
are equal to 1, 3, and 15 for each of them. Figure 1 compares the evolution of all 12 of 
these moments for the Martyna-Klein-Tuckerman model. The moments are fully consistent 
with the ergodic distribution. The moments are reproduced to an accuracy of about four 
significant figures during a simulation of 10^^ timesteps. It should be noted that because 
the last of the MKT equations, ^ — 1 , forces the longtime value of ( ) to be unity, 

the ( numerical ) error of that moment, of order 10“® , is smaller than that for all the rest . 
For all the other data, using a timestep of dt = 0.001 the single-step integration error is of 
order 10“^^ , about the same as double-precision roundoff error. The number of oscillations 
represented is about 10^ , quite consistent with a random error on the order of the inverse 
square root of the number of independent sample oscillations. 
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Second Moments Fourth Moments Sixth Moments 



Time range shown is 100 000 000 to 200 000 000 


Figure 1; Typical variation of cumulative second, fourth, and sixth moments of the MKT oscil¬ 
lator. The exact moments, 1, 3, and 15, are reproduced with four-figure accuracy for simulations 
describing on the order of 10® vibrations. 

B. Chaoticity and the Largest Lyapunov Exponent 


Chaos is an essential ingredient of ergodicity. Chaos can be quantihed by measuring 
the evolution of Lyapunov instability, the ongoing tendency toward the exponentially-fast 
separation of neighboring phase-space trajectories. A steady-state measurement of Lyapunov 
instability can be implemented by forcing a tethered “satellite” trajectory to follow the lead 
of a “reference” trajectory. Both reference and satellite follow exactly the same motion 
equations but with the reference-to-satellite separation continually constrained by rescaling 
its phase-space separation A at the end of every timestep : 


A^ — frify 1 ^t+dt 


A 


t+dt 


lA 


t-\-dt \ 


The rescaling of the separation, per unit time, dehnes the local Lyapunov exponent \{t) : 

\^t+dt\ 


\(t) = {1/dt) In 


lA, 


~ {1/dt) ln[ e 


\dt 1 _ 


= A . 
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Figure 2: Instantaneous values of the Lyapunov exponent ( on the left ) and a histogram of 
integrated million-timestep averages ( on the right ). The Gaussian distribution at the right is 
drawn with the observed mean and standard deviation and is a near-perfect fit to the data. These 
are probability densities so that the vertical scales are set by normalization. 

The long-time-averaged value of this “local” time-dependent Lyapunov exponent is the ex¬ 
ponent A : 

{ X{t) ) = A ~ ^ 111 

The continuous limit dt —)■ 0 can be imposed by using an appropriate Lagrange multiplier—. 

Figure 2 compares a histogram of 10,000 values of the instantaneous “local” exponent 
\(t) , separated by 1000 timesteps with dt = 0.001 , to a histogram of integrated averages. 
Each averaged exponent represents one million timesteps, At = 1000 . The reference-to- 
satellite offset vector A has a length 0.000001 . The time averages have a mean and standard 
deviation : 

{ A(f) )At=iooo = 0.066 ± 0.011 . 

If the distribution were Gaussian the probability of Ending a vanishing integrated exponent 
in a million trials would be about ~ 10“® . The Gaussian shown in the Figure leads 


At+dt 
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to two conclusions: [ 1 ] one million timesteps are clearly sufficient for the Central Limit 
Theorem to apply, so that [ 2 ] the likelihood of hnding a false-negative time average is 
indeed about one in 100 000 000. For all three of the time-reversible harmonic oscillator 
thermostat models, HH, JB, and MKT, samples of one million time-averaged Lyapunov 
exponents were examined^i. The initial conditions were chosen randomly from the four¬ 
dimensional stationary distributions. The data were consistent with chaos and with the 
absence of regular toroidal trajectories ( which would correspond to vanishing Lyapunov 
exponents ). These Lyapunov exponent investigations establish that the measure of any 
nonergodic component is less than 0.000001 . 

C. Analysis Near the Two MKT Fixed Points: (0,0,-l,-|-l) and (0,0,-|-l,-l) 

The MKT oscillator has two separated hxed points where the coordinate and momentum 
vanish, q = p = 0 . The thermostat variables are (C, 0 = (“1) +1) (+1) “1) • The appar¬ 

ent ergodicity of the oscillator implied by the Gaussian moments and the positive Lyapunov 
exponent suggests that neither hxed point is stable. To show that both are actually expo¬ 
nentially unstable we linearize the equations of motion about the hxed points by considering 
a perturbation vector 6 = {6q, Sp, 5^, S^) . 

We begin with the hxed point singled out for analysis by Patra and Bhattacharya^ 

{(1,pX,0 = (0,0,-l,+l) : 

Sq = 5p ] 6p = -6q + 6p ] X = X 4 = “25^ . 

Another time diherentiation provides the separated equations of motion for the perturbations 
in the {q,p) and (C,0 planes. The {q,p) perturbations are linearly unstable, both in the 
same manner : 

5q 6q -|- , 6p 5p -|- 6p . 

The perturbations in the (C,0 plane are stable, again in the same manner : 

X — — 25 ^ — ( 5 ^ ; 4 ~ — 4 • 

Evidently the (C;0 toward this hxed point complements the corresponding Lyapunov- 
unstable exit how in the {q, p) plane. 
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Exactly similar analysis can be carried out at the other fixed point {q,p,C^O — 

(0,0,+1,-1) : 

= ~^q ~ ^p ] ~ ~ +2(5^ . 

This time perturbations {6q,6p) in the {q,p) plane are linearly stable rather than unstable, 
and both in the same manner : 

Sq = -Sq - 6q ] Sp = Sp - 6p . 

In parallel, the perturbations in the (C,0 plane are unstable : 

'S(^ = — 26(^ + = — 2 ( 5 ^ + < 5 ^ . 

In summary, both the fixed points are exponentially unstable, with stable entrance and 
unstable exit flows balancing in the steady state. 

In addition to the mirror symmetry mentioned in the first Section all three sets of ther- 
mostated oscillator equations exhibit a time-reversal symmetry in which the signs of {p, 
and the time all change while that of the coordinate q does not. This symmetry implies that 
the attractors and repellors change roles in the time-reversed dynamics, with damped sta¬ 
ble oscillation reversed to give unstable divergent oscillation and vice versa. In either case 
the oscillator equation 6q = —6q ± 6q corresponds to successive amplitude changes larger 
or smaller by a factor 6.1 . The parallel thermostat equation near the (C^O fixed points, 
(5^ = —2(5^ ± corresponds to successive amplitude changes of a factor of 3.3 at the control 
variable’s turning points, smaller because the characteristic frequency of these thermostat 
variables is greater. 

Evidently the {q,p) flow toward the fixed point (0,0) competes with the exit unstable 
flow in the {(, plane. One might expect that exponential divergence would overwhelm 
the exponential slowing. In order better to understand the fixed point flows we collect 
trajectory points in three thin four-dimensional slabs centered on {q,p) = (0,0) and on 
(C)0 = (“1)+1) aiid (+1,-1) . The thickness of the three slabs are all 10“®/^ . See 
Figure 3 . The two cross-sectional views of {q,p) look precisely similar at the two (C^O 
“fixed points”. While q and p are both small ( corresponding to a measure of 0.00001 ) 
the (C, 0 flow parallels the curve joining the two fixed points and emphasized in the center 
of the Figure. The relatively lengthy (C,0 excursions correspond to much smaller {q,p) 
tracks. In the thin slabs with ~ — 1 , the coordinate changes by more than a factor of six 




Figure 3: {q,p) data are plotted for values close to the “attractor” ( left ) and “repellor” ( right ) 

values of {C,C) =(“1)+1) and (+1,-1). The middle panel shows (C,0 data for {q,p) values close 
to (0,0) . Both “fixed points” are exponentially unstable. The ranges shown for all four variables 
are from -4 to +4. 

between crossings, and the amplitude of the {q,p) motion is much less. These two effects are 
responsible for the misleading appearance of “holes” in the {q,p) projections. We will soon 
show that the density in the full four-dimensional {q,p,CyO space is actually completely 
uniform near both of the hxed points. It is simply the jumps in q coupled with the slow flow 
in p that accounts for the low-density appearance emphasized by Patra and Bhattacharyai^. 

D. Exponential Motion Near the Fixed Points 

Near the two hxed points the how is dominated by the source-to-sink S-curve shown in 
the central panel of Figure 3 and well approximated by the (C,0 projection in the right 
panel of Figure 4. It is educational to conhrm the linear stability analysis by considering 
the how shown in Figure 4, starting very near the (C,0 “source” and {q,p) “sink” : 

{q,p, C, 0 = 10-^2^ +1 + 10-l^ -1 - 10-12). 

At the right in Figure 4 we plot the {q,p) and (C^O trajectories. From a visual standpoint 
the (C, 0 trajectory leaves the repellor at a time just past 40 and quickly moves to the attac- 
tor, settling there near a time 60, and remaining there until just past 140. During all that 
time the {q,p) coordinates appear motionless. Their distance from the origin corresponds 
to the heavy line in the left panel. The distance to the {(, repellor is the light line there. 
The medium left-panel line is the distance to the {(,0 attractor reached near time 60. 

After a time of 148 chaos ensues and the linear stability analysis visible in the left panel 
no longer applies. The oscillator and thermostat plane trajectories in the right panel are 


9 









Figure 4: The left panel shows the {q,p) and (CjO separations from (0,0) and (=Fl,il) for a 
simulation started very near (0,0,+1,-1) . The heavy blue line shows the {q,p) separation, nearly 
invisible until a time of 140, then spiraling away from the origin while the light and medium red 
lines show that the (C, 0 track travels from the repellor ( along the light line ) and toward the 
attractor ( medium ) by a time of 60, remaining near there until visible chaos ensues at time 140 . 

typical and show that the {q,p) motion is slower and less vigorons than the (C,0 motion. 

E. Probability Density Near the Two Fixed Points 

In Figure 5 we plot ( on logarithmic scales ) the nnmber of points ont of 10^° ( lower 
curve ), 10^^ , and 10^^ ( upper curve ) lying within a distance r of the two hxed points, 
with ln(r) ranging from -5 to + 2. Because dln(r) = {dr/r) the density in four-dimensional 
space should vary as rather than , which it does, very accurately. The measure at the 
hxed points is equal to {l/47r‘^e)dqdpd(d^ . Because the volume of a four-dimensional sphere 
of radius r is {7T‘^r‘^/2) the probability of hnding a trajectory point near one of the two hxed 
points within the smallest radius r = e~^ is e“^^/8 ~ 9 x 10“^^ , explaining why such points 
occur rarely even on a 10^^-point trajectory, as is shown in the Figure. 
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Figure 5: The number of trajectory points in the vicinity of the attractor and repellor are shown 
for simulations of 10^*^ , 10^^ , and 10^^ points. The two curves for each simulation are only 
distinguishable for small r where their fluctuations are visible. The data are consistent with a 
Gaussian density and show that the slope near 4 persists up to spherically-averaged r values of 
order unity. 

F. Grid-Based Measures 

We can also confirm the Gaussian solution of the MKT equations by computing the 
measures of 81 four-dimensional nonoverlapping balls of radius 0.50 arranged on a hypercubic 
3 X 3 X 3 X 3 grid centered on the origin. The measures of the balls are inversely proportional 
to the number of nonzero exponents. The ball at the origin has none. There are respectively 
8, 24, 32, and 16 balls with their centers having 1, 2, 3, and 4 nonvanishing exponents. 
A simulation with 10^^ timesteps gave measures of 0.00719, 0.00445, 0.00276, 0.001706, 
0.001056 for the hve ball types. Each of these measures is close to 0.620 times that of its 
predecessor, as expected for product measures of four independent Gaussians. The two of 
the 81 balls centered on the dynamics’ unstable hxed points show nothing out of the ordinary 
in their measures relative to the other 22 two-exponent balls. This statistical test, which 
could be rehned indehnitely in complexity, is, like all the rest, consistent with ergodicity of 
the MKT oscillator equations. 
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III. SUMMARY 


The thermostated oscillator model introduced by Martyna, Klein, and Tuckerman, and 
explored here in more detail, is one of three simple systems exhibiting a smooth Gaussian 
distribution accompanied by a complex chaotic dynamics. All three are important from the 
pedagogical standpoint, and are also useful as thermostats generating all of Gibbs’ canonical 
distribution. 

Here we have summarized the ( compelling ) evidence for the ergodicity of the MKT 
oscillator in order to close out the competition for the 2014 Snook Prize. We have used six 
different and independent methods to assess the ergodicity of the MKT oscillator: [1] the 
moments of the Gaussian distribution; [2] the chaos, as opposed to regularity, of billions 
of independent trajectories; [3] the instability of the flows near both hxed points; [4] the 
exponentially growing separation from both fixed points; [5] the uniform probability density 
in the vicinity of these unstable hxed points and [6] the expected relative measures within a 
set of 81 hyperspheres centered on the lattice nodes of a four-dimensional hypercubic lattice. 

All of these methods reach the same conclusion, that solutions of the coupled equations 
are ergodic. We hope that this summary article will prove useful to investigators of ergodicity 
in other simple dynamical systems. 

In view of the very intricate Lyapunov instability of the Martyna-Klein-Tuckerman sys¬ 
tem this ergodic Gaussian distribution is outstanding in its simplicity. In view of the con¬ 
tributions of Puneet Patra and Glint Sprott to the understanding of this problem we have 
divided the 2014 Snook Prize equally among ourselves and themselves. We intend to for¬ 
mulate another Snook Prize problem in the summer of 2015 and would be very grateful for 
suggestions from the readers. We thank Puneet Kumar Patra and Julien Glinton Sprott for 
helpful support and Ben Leimkuhler and Mark Tuckerman for stimulating comments. 
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